Effects of Discrete Breathers on Heat Conduction 
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Intensive studies in the past decades have suggested that the heat conductivity k diverges with the 
system size L as k ~ L a in one dimensional momentum conserving nonlinear lattices and the value 
of a is universal. But in the Fermi-Pasta-Ulam-,9 lattices with next-nearest-neighbor interactions 
we find that a strongly depends on 7, the ratio of the next-nearest-neighbor coupling to the nearest- 
neighbor coupling. We relate the 7-dependent heat conduction to the interactions between the 
long-wavelength phonons and the randomly distributed discrete breathers. Our results provide an 
evidence to show that the nonlinear excitations affect the heat transport. 
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In the studies of non-electronic heat conduction, it is 
an important progress to realize that the heat conduc- 
tivity k diverges with the system size L as k ~ L a in 
one dimensional momentum conserving nonlinear chains 
[l|, [2j. The value of the exponent a is believed to be 
universal [3j|, though there are lots of debates on what 
value(s) it takes, (e.g., if there exists one universal class 
with a = I d H or two with a = \ and a = § [l-[i| 
is still controversial.) The universality of a roots in the 
theory first proposed by Peierls [9|, where the essence of 
the non-electronic heat conduction at low temperatures 
is modeled as a weakly interacting phonon gas. Based on 
this model, a universal heat conduction law constrained 
only by the dimensionality of a system regardless of the 
details of its microscopic dynamics is thereby expected. 

On the other hand, at high temperatures nonlinear ex- 
citations such as traveling solitary waves [10| and discrete 
breathers (DBs) [11| are ubiquitous in nonlinear lattice 
systems. Hence interactions between phonons and non- 
linear excitations should be studied and taken into ac- 
count when their effects are considerable. As nonlinear 
excitations involves microscopic dynamical details, if a 
universal heat conduction law still exists certainly de- 
serves careful investigations. In this respect quite few 
studies have been reported. Early work of our group 
showed that traveling solitary waves may play an impor- 
tant role in heat conduction of the Fermi-Pasta-Ulam- 
/3 (FPU-/3) chain s 112. 1 131 . but this was argued against 



by some authors 



1, 



m- 



Besides, DBs have also been 



proposed as a phonon scattering mechanism |16| for the 
normal heat conduction numerically observed in the har- 
monic chains with nonlinear on-site potentials |17| and 
in the rotator chains [18|. In spite of these studies, at 
present whether — and if yes how the nonlinear excita- 
tions would affect the heat transport in low dimensional 
momentum conserving systems is still an open question. 

In this Letter we present a study to show how DBs play 
a role in the heat conduction in one dimensional nonlin- 
ear lattices. Specifically, we perform numerical analysis 



to investigate the heat conduction in FPU-/3 chains with 
the next-nearest-neighbor (NNN) coupling. Our results 
show clearly that a varies continuously with the ratio of 
the NNN coupling to the nearest-neighbor (NN) coupling, 
suggesting our system does not belong to any universal- 
ity class characterized by a constant a. Moreover, we 
find a is correlated to the overlap of the phonons' spec- 
tra and the DBs' spectra, and in such a sense we relate 
the 7-dependent heat conduction behavior to the phonon 
scattering by the DBs. 

Our system is a chain of N identical particles with both 
the NN and NNN interactions [19] whose Hamiltonian is 

2 
H = Efi + V (««+i - *) + 1 V («+a - %)]• (1) 

i 

Here qi is the displacement of the ith particle from its 
equilibrium position and pt its momentum, and the po- 
tential is of the FPU-/3 type; i.e., V (x) = \x 2 + \x 4 . 
Both the mass fi and the lattice constant are set to be 
unit. The parameter 7 is tunable; it specifies the compar- 
ative strength of the NNN coupling. 7 = corresponds 
to the conventional FPU-/? system. 

We employ the reverse non-equilibrium molecular dy- 
namics simulation method (RNEMD) [20( to build the 
nonequilibrium stationary state across the system. Com- 
pared with the usual method that couples the system 
to two heat baths at its two ends, the RNEMD is 
advantageous in suppressing the boundary effects and 
therefore leads to a faster convergence to the stationary 
state. Meanwhile the RNEMD keeps the total energy 
and momentum of the system unchanged. Following the 
RNEMD [20||, first the periodic boundary condition is 
imposed; then the chain (now a circle) is divided into M 
slabs each contains n — -tt particles. The instantaneous 
local kinetic temperature T/. in slab k (k = 1, • • • , M) 
is defined by T k = -^ Y%t n {k-i)+iPh where k B is the 
Boltzmann constant (set to be unit) and the sum goes 
over all n particles in slab k. The temperature profile of 
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Figure 1: (Color online) Top panels: the temperature profile 
for 7 = (a) and 7 = f (b) with the effective system size 
L = 2496. Bottom panels: the heat conductivity n versus L 
for 7 = (c) and 7 = 1 (d). The dashed lines are for the best 
fitting of k ~ L a , suggesting a = 0.325 ± 0.002 for 7 = (c) 
and a = 0.35 ± 0.02 for 7 = 1 (d). 



the stationary state can be represented by the time aver- 
age of T fc (fc = 1, ■ • • , M), denoted by (T k ). Next, slab 1 
is set to be the cold slab and slab -^ + 1 as the hot slab. 
The heat flux is generated by exchanging the momenta 
of the hottest (coldest) particle in the cold (hot) slab at a 
fixed frequency / cxc (the exchange frequency) . This pro- 
cedure results in a redistribution of a certain amount of 
kinetic energy AE — ^2 h(Ph ~ Pc) during time t. Here 
the subscript h and c refer to the hot and cold particles 
whose momenta are interchanged, and the sum takes all 
exchanges in time t. In the stationary state the relax- 
ation of AE will drive two heat fluxes of J — ^- across 
the system, because due to the periodic boundary condi- 
tion we have in fact two identical chains between the hot 
and cold slabs. L = y — n is the effective length of the 
two chains (deducting the hot and cold slabs) and J is 
the heat flux crossing each of them. Once the stationary 
state is reached, a temperature gradient between the hot 
and cold slabs is expected and the thermal conductivity 
k can then be measured by assuming the Fourier law; i.e., 
- - -I_ 

n, — VT . 

We start our simulations with a fully thermalized chain 
at temperature T = 2.5. The velocity- Ver let algorithm 
[21| with a time step 0.01 is used to evolve the system, 
and M = 80 and / cxc = 0.1 are adopted for the RNEMD. 
For each system size a transient stage of time 10 6 is dis- 
carded (which has been verified to be long enough for 
reaching the stationary state); then the next evolution 
of time 10 7 is performed for the time average. We have 
checked that our results do not qualitatively depend on 
the particular parameter values taken here. 

Before presenting our main results, it is interesting to 
make a quick comparison between our simulations and 
those by different methods for 7 = [5j and 7 = 1 [19J . 
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Figure 2: (Color online) The dependence of a on the param- 
eter 7. The vertical dashed line indicates 7 C = 0.25. Error 
bars give the standard error for evaluating a by linearly fixing 
log k versus log L. 



For 7 = 0, i.e. the conventional FPU-/3 system, Fig. 1(a) 
shows the temperature profile for L = 2496; it can be seen 
that a constant temperature gradient is well established 
between the hot and cold slabs. In addition, for larger 
system sizes the temperature profiles (not shown) have 
been checked to be the same upon a rescaling. Fig. 1(c) 
shows the dependence of k on the effective system size L; 
it suggests that k diverges as ~ L a with a = 0.325±0.002 
which we emphasize to be very close to the predicted 
value I by the hydrodynamic theory_|4| and the result 
of a recent careful numerical study [5|. For 7 = 1 we 
have obtained the similar results as 7 = [see Fig. (b) 
and (d)], but the best fitting performed over 2496 < L < 
19968 gives a = 0.35 ±0.02 instead [see Fig. 1(d)]. Note 
that our a value for 7 = 1 is remarkably different from 
that given in Ref. [19j where a was evaluated over much 
shorter system sizes (L < 2000) and therefore fails to 
capture the divergence of k in the thermodynamical limit. 

The first main result of our study is shown in Fig. 2 
for the divergence exponent a versus the parameter 7. It 
can be seen that as 7 changes from to 1, a decreases 
and reaches its minimum a m i n s* 0.24 at 7 w r y c = 0.25, 
then increases up to about 0.35 at 7 = 1 with a trend of 
saturation [22]. The fact that a changes continuously in 
a range is in clear contrast to the existence of the general 
a value(s) independently of the dynamics. 

As phonons are the heat energy carriers in our lat- 
tices, a takes its minimum at j c implies that around 7 C 
the NNN coupling may enhance the phonon scattering 
and thus give rise to smaller values of a. To clarify this 
point the phonon dispersion relation turns out to be very 
suggestive. Keeping the harmonic coupling in both the 
NN and NNN interactions, the dispersion relation reads 
uj q = 2 [sin I + 7 sin q]? , where q is the wave number 
and ujq the corresponding frequency. Interestingly, j c is 
a transition value for the phonon dispersion relation as 
well (see Fig. 3): For 7 < j c the maximum frequency 
cjjr = 2 corresponds to the boundary of the Brillouin 




Figure 3: (Color online) Phonon dispersion relation for the 
FPU-/3 system with the NNN interactions. The curves from 
bottom to top correspond to 7 = 0, 0.25, and 1 respectively. 
The horizontal dash-dotted line indicates uj n = 2. 



zone q = n but for 7 > j c it grows larger than w, and 
moves away from the boundary. For 7 = j c the group 
velocity is close to zero in a wider q region. This property 
favors the DBs in the presence of the nonlinearities [23| 
which we will discuss later [see Fig. 5(b)]. 

Since the work by Peyrard 24| it has been known that 
the temperature activated DBs may be crucial for the 
energy transport and other dynamical processes. Inter- 
esting examples include the meltin g tr ansitions in solids 
and folding in polypeptide chains |25| . To study if the 
DBs may have any effects on the heat conduction in our 
system, it is necessary to check if DBs exist at the focused 
temperature T — 2.5. We apply the method in Ref. [17|; 
i.e., a chain of TV = 2000 particles is first thermalized at 
temperature T — 2.5, then the heat baths are removed 
and the absorbing boundary conditions are imposed to 
the two ends of the chain [26j . After all the mobile exci- 
tations such as phonons and solitary waves are absorbed, 
a few standing breathers may emerge in the internal seg- 
ment of the chain. The snapshot of the energy profile 
after a long time (8 x 10 5 ) absorbtion is presented in Fig. 
4(a) and (b) for 7 = and 7 = j c , respectively; in both 
of them the DBs can be well recognized. We have verified 
that this is also the case for other 7 values in [0, 1]. 

Now we study the interactions between the DBs and 
the phonons. For this aim we calculate the power spectra 
P(oj) of the residual thermal fluctuations after long time 
absorbing. To facilitate the computation short chains of 
size N — 200 are considered. The results for 7 = and 
7 = 7 C are plotted in Fig. 4(c) and Fig. 4(d) for a com- 
parison; It shows that for 7 = the DBs' frequencies 
are outside the linear phonon band of < u> < u>^, in 
agreement with the classical DB theory 2JJ • Due to this 
frequency mismatch we may conjecture the lack of inter- 
actions between the DBs and the phonons, which in turn 
implies the null effects of the DBs on the heat conduc- 
tion. However, in clear contrast, for 7 = j c a portion 
of the DBs' frequencies appear inside the linear phonon 
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Figure 4: (a) and (c): The snapshot of the energy distribution 
and the power spectrum of the residual thermal fluctuations 
for 7 = 0. (b) and (d): the corresponding results for 7 = 0.25. 
While for 7 = no in-band components can be identified, 
there are considerable in-band components for 7 = 0.25 [see 
the insert for a zoom of the boxed in-band region in (d)]. 



band, suggesting the existence of the in-band DBs [25 
|30J. This is well shown in the insert of Fig. 4(d) where 
the collective modes in the linear phonon band can be 
clearly recognized. As the in-band DBs can randomly 
distribute along the lattice and interact with phonons, 
they introduce an inherent disorder [27] and play roles of 
random scatters to the phonons. This may explain why 
a is smaller in the case of 7 = j c . If this is true, combin- 
ing the results given in Fig. 2 we may expect that as 7 is 
increased, the interactions between the in-band DBs and 
the phonons would become stronger and stronger (weaker 
and weaker) for < 7 < j c (^ c < 7 < 1). To check if this 
is the case, we define e = L " P(ui)dui/ L P(u>)du>, the 
ratio of the energy of the collective modes within the lin- 
ear phonon band to the total energy of the residual ther- 
mal fluctuations, as a measure of the interaction intensity 
between the in-band DBs and the phonons. For several 
typical 7 values we first calculate the power spectra in the 
same way as in Fig. 4(c) and (d), then evaluate e and 
summarize the results in Fig. 5(a). It shows that e ver- 
sus 7 is in good consistence with a versus 7 (see Fig. 2). 
This consistence verifies that the heat conduction behav- 
ior is indeed dependent on the DB-phonon interactions 
in our system. (Note that in Fig. 5(a) the maximum of e 
does not correspond to 7 C exactly but a slightly smaller 7 
value; this discrepancy may be a result of the big statis- 
tical errors in the evaluated power spectra of the residual 
thermal fluctuations where short chains of TV = 200 have 
to be used on account of computation cost.) 

In the following let us turn to the question why the 
DB-phonon interactions are the strongest at 7 = 7 C . Our 
study [see Fig. 5(b) later] suggests the DBs' concentra- 
tion is the highest at 7 = j c , which we conjecture results 
in the highest probability for the DBs to interact with 
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Figure 5: (Color online) (a) e versus 7 and (b) ft versus 7 for 
e(0) = 0.005 (squares), 0.045 (dots) and 0.18 (triangles). The 
vertical dashed line indicates 7 C = 0.25 in both panels. 



the phonons. The concentration of DBs is approximately 
e k B T in an equilibrium state, where e s h is the energy 
threshold for creating the DBs [27]. It implies that e s h is 
the smallest if the DBs' concentration is the highest (at 
7 = 7 C ). To check if this is the case we consider an energy 
relaxation process. Initially we give a kinetic energy e(0) 
to the particle centered at the chain of TV = 2000, a zero 
velocity to all others and a zero displacement to all parti- 
cles. Then absorbing boundary conditions are imposed, 
and after a transient time of 2 x 10 5 we calculate e(t), 
the total energy remained in the chain at time t, up to 
t = 8 x 10 5 . We find it decays as e{t) ~ t~ p J23| and the 
exponent ft depends on both the initial excitation energy 
e(0) and 7: When e(0) is small enough harmonic behav- 
iors manifest themselves and ft — > 0.5 as being pointed 
out in [31]; but when e(0) is large enough long-lived DBs 
could form and thereby ft — > 0. We increase e(0) progres- 
sively and find it is exactly at 7 = j c where the signal 
for the DBs appear first [see Fig. 5(b)]. Hence we may 
conclude for 7 = j c the value of e s h is the smallest and 
the concentration of the DBs is the highest. 

In summary, we have studied a one dimensional lat- 
tice of the FPU-/3 type with the NNN coupling. We have 
shown that tuning the NNN coupling may change the 
concentration of the DBs, which in turn vary the DB- 
phonon interaction intensity, and consequently affect the 
heat conduction characteristics of the system. Depend- 
ing on the NNN coupling, the divergence exponent a of 
the heat conductivity may continuously change from 0.24 
to about |. In contrast to the previously suggested gen- 
erality class (es), our study reveals a new regime in one 
dimensional heat conduction featuring the interactions 
between nonlinear excitations and phonons. At present 
the understanding to this new regime is still primitive; 
further investigations are necessary and desired. 
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